function [psi,rhoj,sj,sHj, sNHj, s0j,s0Hj, s0NHj, etaEj,etaDj,piDprimej,dealer_idsSQ,uB,uE,piB,piE,YB,YE,Ecost,  q_data, rhoj_data, s0j_data,xij_data, YBno, EUno, uE2, psi_opt, rhoj_opt, vdjSQout, W, Wxij, Wvdj, WxijH, W_shrink, WwH, WxijwH, WvdjwH, WxijHwH, W_shrinkwH] = status_quo_fun(SQ, CF, side, access, Yg, new_periodg, new_period, vdjSQ, estimate, Mhat)
   
global T Imax Dmax xij_for_investor homedealer eta loyalty

if estimate ==1
    K=size(new_periodg,1); 
else 
    K=T;
end 

kc=27:29;

%Initialize  
psi       = NaN* ones(K,Imax);
psi_opt   = NaN* ones(K,Imax);
rhoj      = NaN* ones(K,Imax);
rhoj_opt  = NaN* ones(K,Imax);
sj        = NaN* ones(K,Imax);
sHj       = NaN* ones(K,Imax);
sNHj      = NaN* ones(K,Imax, Imax);
s0Hj      = NaN* ones(K,Imax);
s0NHj     = NaN* ones(K,Imax, Imax);
s0j       = NaN* ones(K,Imax);
etaEj     = NaN* ones(K,Imax);
etaDj     = NaN* ones(K,Imax);
piDprimej = NaN* ones(K,Imax);
dealer_idsSQ= NaN* ones(K,Imax);
uB        = NaN* ones(K,Imax);
uE        = NaN* ones(K,Imax);
uE2       = NaN* ones(K,Imax);
piB       = NaN* ones(K,Imax);
piE       = NaN* ones(K,Imax);
YB        = NaN* ones(K,Imax);
YE        = NaN* ones(K,Imax);
Ecost     = NaN* ones(K,Imax);
YBno      = NaN* ones(K,Imax);
EUno      = NaN* ones(K,Imax);
vdjSQout  = NaN* ones(K, Imax);
W         = NaN* ones(K, 1);
W_shrink  = NaN* ones(K, 1);
Wxij      = NaN* ones(K, 1);
WxijH      = NaN* ones(K, 1);
Wvdj      = NaN* ones(K, 1);

WwH         = NaN* ones(K, 1);
W_shrinkwH  = NaN* ones(K, 1);
WxijwH      = NaN* ones(K, 1);
WxijHwH      = NaN* ones(K, 1);
WvdjwH      = NaN* ones(K, 1);

q_data     = NaN* ones(K,Imax);
rhoj_data  = NaN* ones(K,Imax);
s0j_data   = NaN* ones(K,Imax);
xij_data   = NaN* ones(K,Imax);


for k=1:K

    %find indices
    date    = new_period(k,1);     
    t       = find(new_periodg(:,1)==date);

    %If period non-missing and (below) parameters non-missing
    if ~isempty(t)
    idxt    = [new_periodg(t,2):new_periodg(t,3)]'; %lines in Y
        
    %label
    Yt      = Yg(idxt,:);
    
    %If retail investors get no home dealer benefit
    %if side==1
    %    Yt(Yt(:,1)==2,24)=Yt(Yt(:,1)==2,25);
    %else
    %    Yt(Yt(:,1)==4,24)=Yt(Yt(:,1)==4,25);
    %end

    if estimate==0
        paramst = Yg(new_periodg(t,2),kc);
    else
        paramst=Mhat(t,:);
    end 
        

    if ~isnan(sum(paramst)) 
   
    %find the line with each dealer's first trade
     l=1;
     for j=1:length(Dmax)
        if sum(Yt(:,16)==Dmax(j))>0
            idxj(l)=find(Yt(:,16)==Dmax(j), 1, 'first');
            l=l+1;
        end
     end
    
    q0      = Yt(idxj, 18);
    rhoj0   = Yt(idxj, 13);
    
    if estimate==0
        vdjSQt  = vdjSQ(t,:);     
    else
        vdjSQt  = []; 
    end 
        
    if CF<0, q0=vdjSQt';
    end 

    if  loyalty==0 %set home dealer benefit to zero
        Yt(:,24) =Yt(:,23);
    end 
    
    s0j0    = Yt(idxj, 12);
    xij0    = Yt(idxj, 23); %basequality 


    if estimate==1
       
        %EXTENSION WITH HOME DEALER
        if homedealer==1 && eta==0
            [OFS, psiSQt, rhojSQt, sjSQt, sHjSQt, sNHjSQt,  s0jSQt, s0HjSQt, s0NHjSQt, etaEjSQt, etaDjSQt, vdjSQtout,piDprimejSQt, ...
            dealer_idsSQt, uB_SQt, uE_SQt,piB_SQt,piE_SQt,YB_SQt, YE_SQt, Ecost_SQt, uE2_SQt, psiSQ_optt, rhojSQ_optt, Wt, Wxijt, Wvdjt,WxijHt,  W_shrinkt, WtwH, WxijtwH, WvdjtwH,WxijHtwH,  W_shrinktwH ] = model_prediction_avg_H(q0, Yt, paramst, side, SQ, CF,access, [], idxj);                                                                                          
            
        %BENCHMARK MODEL:
        elseif homedealer==0 && eta==0
            [OFS, psiSQt, rhojSQt, sjSQt, s0jSQt, etaEjSQt, etaDjSQt, vdjSQtout,piDprimejSQt, ...
            dealer_idsSQt, uB_SQt, uE_SQt,piB_SQt,piE_SQt,YB_SQt, YE_SQt, Ecost_SQt, uE2_SQt, psiSQ_optt, rhojSQ_optt, Wt, Wxijt, Wvdjt, W_shrinkt] = model_prediction_avg(q0, Yt, paramst, side, SQ, CF,access, [],idxj, t);                                                                                          
       
        %ROBUSTNESS: Eta>0
        elseif homedealer==0 && eta>0
            [OFS, psiSQt, rhojSQt, sjSQt, s0jSQt, etaEjSQt, etaDjSQt, vdjSQtout,piDprimejSQt, ...
            dealer_idsSQt, uB_SQt, uE_SQt,piB_SQt,piE_SQt,YB_SQt, YE_SQt, Ecost_SQt, uE2_SQt, psiSQ_optt, rhojSQ_optt, Wt, Wxijt, Wvdjt, W_shrinkt] = model_prediction_avg_eta(q0, Yt, paramst, side, estimate, access, [],idxj, CF);                                                                                          
        end 
            
    else
          
        %EXTENSION WITH HOME DEALER
        if homedealer==1
            
            [OFS, psiSQt, rhojSQt, sjSQt, sHjSQt, sNHjSQt, s0jSQt, s0HjSQt, s0NHjSQt,etaEjSQt, etaDjSQt, vdjSQtout,piDprimejSQt, ...
            dealer_idsSQt, uB_SQt, uE_SQt,piB_SQt,piE_SQt,YB_SQt, YE_SQt, Ecost_SQt, uE2_SQt, psiSQ_optt, rhojSQ_optt, Wt, Wxijt, Wvdjt, WxijHt, W_shrinkt,  WtwH, WxijtwH, WvdjtwH,WxijHtwH,  W_shrinktwH] = model_prediction_avg_H(q0, Yt, paramst, side, SQ, CF,access, vdjSQt, idxj);                                                                                              
             
        %BENCHMARK MODEL: 
        elseif homedealer==0
            [OFS, psiSQt, rhojSQt, sjSQt, s0jSQt, etaEjSQt, etaDjSQt, vdjSQtout,piDprimejSQt, ...
            dealer_idsSQt, uB_SQt, uE_SQt,piB_SQt,piE_SQt,YB_SQt, YE_SQt, Ecost_SQt, uE2_SQt, psiSQ_optt, rhojSQ_optt, Wt, Wxijt, Wvdjt, W_shrinkt] = model_prediction_avg(q0, Yt, paramst, side, SQ, CF,access, vdjSQt,idxj, t);                                                                                          
        
        %ROBUSTNESS: Eta>0
        elseif homedealer==0 && eta>0
            [OFS, psiSQt, rhojSQt, sjSQt, s0jSQt, etaEjSQt, etaDjSQt, vdjSQtout,piDprimejSQt, ...
            dealer_idsSQt, uB_SQt, uE_SQt,piB_SQt,piE_SQt,YB_SQt, YE_SQt, Ecost_SQt, uE2_SQt, psiSQ_optt, rhojSQ_optt, Wt, Wxijt, Wvdjt, W_shrinkt] = model_prediction_avg_eta(q0, Yt, paramst, side, estimate,access, vdjSQt,idxj, CF);   
        end 

    end
       

    %YIELDS
     idxd       = find(ismember(Dmax,dealer_idsSQt));
     YBnot      = NaN*ones(Imax,1);
     EUnot      = NaN*ones(Imax,1);
          
      if side==1 
            YBnot(idxd) = Yt(find(Yt(:,17)==1), 5) + paramst(:,1) - xij0.*( (xij_for_investor==1 & homedealer==0) | (homedealer==1 & xij_for_investor==0) );  %theta + mu - basequality
      elseif side==-1 
       YBnot(idxd) = Yt(find(Yt(:,17)==1), 5) + paramst(:,1) + xij0.*( (xij_for_investor==1 & homedealer==0) | (homedealer==1 & xij_for_investor==0) );  %theta + mu + basequality
          
     end 
     EUnot(idxd) = xij0.*( (xij_for_investor==1 & homedealer==0) | (homedealer==1 & xij_for_investor==0) );
     

    %STACK AND STORE 
    idxd            = find(ismember(Dmax,dealer_idsSQt));
    psi(k,:)        = psiSQt'; 
    rhoj(k,:)       = rhojSQt';
    rhoj_opt(k,:)   = rhojSQ_optt';
    psi_opt(k,:)    = psiSQ_optt'; 
    sj(k,:)         = sjSQt';
    if homedealer==1
        sHj(k,:)    = sHjSQt';
        sNHj(k,:,:) = sNHjSQt';
        s0Hj(k,:)   = s0HjSQt';
        s0NHj(k,:,:)= s0NHjSQt';
    end 

    s0j(k,:)        = s0jSQt';
    etaEj(k,:)      = etaEjSQt';
    etaDj(k,:)      = etaDjSQt';
    piDprimej(k,:)  = piDprimejSQt';
    dealer_idsSQ(k,:) = dealer_idsSQt';
    uB(k,:)         = uB_SQt';
    uE(k,:)         = uE_SQt';
    uE2(k,:)        = uE2_SQt';
    piB(k,:)        = piB_SQt';
    piE(k,:)        = piE_SQt';
    YB(k,:)         = YB_SQt';
    YE(k,:)         = YE_SQt';
    Ecost(k,:)      = Ecost_SQt';
    YBno(k,:)       = YBnot';
    EUno(k,:)       = EUnot';
    vdjSQout(k,:)   = vdjSQtout';
    W(k,:)          = Wt';
    W_shrink(k,:)   = W_shrinkt';
    Wxij(k,:)       = Wxijt';
    Wvdj(k,:)       = Wvdjt';

    if homedealer==1
        WxijH(k,:) = WxijHt';
        WwH(k,:)          = WtwH';
        W_shrinkwH(k,:)   = W_shrinktwH';
        WxijwH(k,:)       = WxijtwH';
        WvdjwH(k,:)       = WvdjtwH';
        WxijHwH(k,:)       = WxijHtwH';
    
    end 
       
    %Put actual data in the same matrix format, create matrices: IxK (with 0's if the dealer was not active)
    q0_fullt                = NaN*ones(Imax,1);
    rhoj0_data_fullt        = NaN*ones(Imax,1);
    s0j0_data_fullt         = NaN*ones(Imax,1);
    xij_data_fullt          = NaN*ones(Imax,1);
    rhoj0_data_fullt(idxd,:)  = rhoj0;
    s0j0_data_fullt(idxd,:)   = s0j0;
    xij_data_fullt(idxd,:)    = xij0;
  
    if CF>=0
        q0_fullt(idxd,:)          = q0;
    end 

    q_data(k,:)     = q0_fullt';
    rhoj_data(k,:)  = rhoj0_data_fullt';
    s0j_data(k,:)   = s0j0_data_fullt';
    xij_data(k,:)   = xij_data_fullt';

    end 
    clear idxj 
    end 
end 

